Numerical Study of the Correspondence Between the Dissipative and Fixed Energy 

Abelian Sandpile Models 
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We consider the Abelian sandpile model (ASM) on the large square lattice with a single dissipative 
site (sink). Particles are added by one per unit time at random sites and the resulting density of 
particles is calculated as a function of time. We observe different scenarios of evolution depending 
on the value of initial uniform density (height) ho = 0, 1, 2, 3. During the first stage of the evolution, 
the density of particles increases linearly. Reaching a critical density p c (ho), the system changes 
its behavior sharply and relaxes exponentially to the stationary state of the ASM with p s = 25/8. 
We found numerically that p c (0) = p 3 and p c (ho > 0) 7^ p s - Our observations suggest that the 
equality p c — p 3 holds for more general initial conditions with non-positive heights. In parallel with 
the ASM, we consider the conservative fixed-energy Abelian sandpile model (FES). The extensive 
Monte-Carlo simulations for ho = 0, 1,2,3 have confirmed that in the limit of large lattices p c {ho) 
coincides with the threshold density p t h(ho) of FES. Therefore, pth{ho) can be identified with p s if 
the FES starts its evolution with non-positive uniform height ho < 0. 
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I. INTRODUCTION 

A long-standing discussion of the comparative criti- 
cal properties of the dissipative Abelian sandpile model 
(ASM) l] and the conservative fixed-energy sandpile 
(FES) [2| has gained recently renewed impetus due to 
the works by Fey et al |i3| ■ The point of discussion in a 
laconic form can be reduced to the single question: given 
the same lattice with open and closed boundary condi- 
tions, whether the stationary density of the dissipative 
ASM p s coincides with the threshold density of the FES 
Pth ? Using large-scale simulations on the square lattice, 
the authors of works 0, Q gave a negative answer to this 
question and supported their numerical findings by exact 
solutions for some graphs of higher symmetry. A more 
detailed answer contains a description of the average den- 
sity of grains p(r) for the given density of added particles 
r of the dissipative sandpile: 



P(r) 



^grains at time t t 

W 2 ' r= iV2' 



(1) 



where N is the linear size of the lattice. It was shown 
in that p(r) for the ASM exhibits a transition at the 
point p c which coincides with the threshold density of 
the FES: r = p c = p t h- For r > p c , the system continues 
evolution and its density tends to the stationary value p s 
of the ASM when r — ► 00. 

In this article, we continue the analysis of the corre- 
spondence between critical points of the ASM and FES. 
Confirming the main result of [3] on the transition point 
at p c — pth, we present arguments against the general 
claim pth 7^ p s for the two-dimensional square lattice. 
The focus of the work will be the dependence of the 
transition point on the initial conditions. Specifically, 



we study the nearly closed sandpile, namely the ASM on 
the square lattice with periodic boundary conditions and 
a single dissipative site. One of the results suggested by 
this study is that there exist initial conditions which lead 
during the evolution to the transition point p c coinciding 
with p s = 25/8, and the subsequent evolution does not 
change this density. On the other hand, our simulations 
of the ASM and FES show that the transition point of 
the nearly closed ASM coincides for large lattices with 
the threshold density pth of the FES for all considered 
initial conditions. To fix our notations, we will specify 
the two-dimensional ASM as in 0], i.e. all stable sites 
have heights hi < 4 whereas sites with hi > 5 are unsta- 
ble and topple. 

Traditionally, empty sites (hi = 0) and negative 
heights are not considered neither analytically nor nu- 
merically because they exist in transient states only and 
vanish when a system approaches criticality. Neverthe- 
less, these sites may affect the evolution of the ASM from 
an initial state to the transition point. The dynamics of 
the nearly closed ASM is very close to that of the FES as 
dissipation through the single site is strongly restricted. 
Non-ergodicity of the FES causing a dependence of the 
threshold density upon the initial conditions has been 
discussed recently in [l2j]. A conclusion from this dis- 
cussion was: "the self-organized criticality exhibited by 
the ASM has nothing to do with phase transitions in the 
FES" [T2|. Being reasonable, this statement does not ex- 
clude however specific initial conditions whose evolution 
leads the FES exactly to the stationary density of the 
ASM Ps = 25/8. 

In what follows, we analyze five initial conditions with 
uniform heights ho — 0, 1,2,3 for all lattice sites beside 
the sink. The density p as a function of the density of 
added particles t behaves differently for different ho and 
has a characteristic kink at p c (ho), also depending on 
ho. For ho = 1, the value p c (l) = 3.12528 ... is close to 
the threshold density of the FES obtained by Fey et al 
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Q. For ho = 0, the density p increases linearly with t, 
reaches the value p s ~ 25/8 and remains constant later. 
In Section 2, we give the definition of the model and re- 
mind the origin of the conjectural exact value p s = 25/8. 
In Section 3 we present details of the numerical procedure 
and describe five different scenarios of evolution for the 
initial conditions ho = 0, 1, 2, 3. The special case ho — 
is considered in more details. 

In parallel with the ASM, we present simulations of 
the FES for finite N x N square lattice with various N 
and demonstrate the critical behavior of the threshold 
density pth{ho)- Namely, we find numerically the finite 
size correction 1/N 1 with critical index 7 = 7 (fro) that 
depends on the initial conditions. 



II. THE MODEL AND CRITICAL DENSITY 

We consider the standard Abelian sandpile model [l[ 
on the N x N square lattice. A peculiarity of our con- 
sideration is the nearly closed boundary conditions: the 
lattice is a torus with a single selected site i$. The height 
Zi at any site beside io takes values 0, 1,2,3,4 in stable 
configurations. Particles are added by one at a random 
site increasing the height at that site by one. If the height 
exceeds 4, then the site becomes unstable and topples, 
transferring one particle to each of four neighboring sites. 
The site io serves as a sink where particles disappear. 
During a long evolution, the system gets eventually into 
a set of recurrent configurations. The theory of recurrent 
states of the ASM was developed by Dhar [7[. One of 
important consequences of this theory is the existence of 
a mapping of the set of recurrent configurations onto the 
set of spanning trees of the same lattice. The mapping 
allows one to find the height probabilities in the recurrent 
state and therefore to find the stationary density p a . Ma- 
jumdar and Dhar @ calculated the probability of height 
h = 1 in the thermodynamic limit 



Pi = - -T. 



(2) 



The probabilities P2 , P3 and P4 = 1 — Pi — P2 — P3 were 
found in 9]: 
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Here /„, v — 1,2 are integrals 

i sin Pi det M v daida2dj3\df32 



I, 



1 



(2tt) 4 Mo D(ax,Pi)D{a2,fa)D(ai+a 3> Px+fay 

(6) 



where 



and Mi , M2 are the matrices 
/ 1 1 

Mi = 



1 

3 e i(/9i+/3a) e i(Q2-&) 

4/tt - 1 e i(«i+«2) 1 

\4/lT- 1 e -'(«i+"2) e 2i«2 
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Later on, an exact relation between P2 and P3 was 
proved in [5j, thereby eliminating one of the two integrals. 
A conjecture on the value of the remaining one, based on 
its numerical evaluation to twelve decimal places, then 
led to the following values for the Pj Q, 



P 2 = t - 



1 1 3 12 



P 3 = 



4 27T 7T 2 7T 3 ' 

3 1 12 



8 ' 7T 7T 3 ' 

and for the stationary density, 

Ps = Pi+ 2P 2 + 3P 3 + 4F 4 = 



(10) 
(11) 

(12) 



D(a, /3) = 2- cos(a) - cos(/3) 



(7) 



This value of p s was first suggested by Grassberger [ll| 
based on a high precision calculation of the integrals 
Ii,l2- The present numerical accuracy reaches 10 -25 , 
and forms the basis for the conjecture that p s = 25/8 
exactly. 

Another conjecture, coined in Q as the "density con- 
jecture", identifies p s with the threshold density of the 
FES, pth 0- The threshold density is defined as the 
average maximal number of randomly dropped particles 
per site which allows the relaxation of the closed sand- 
pile to stop. The value pth marks the border of stabil- 
ity: for any p > p t h, an avalanche process started by 
adding a particle never stops with probability one. The 
conjecture pth = Ps was supported by numerical calcula- 
tions [2]. However, recently, Fey et al 0] have obtained 
Pth = Ps + A, where A = 0.000288 to six decimal places. 
Together with the claim by Park [TJ| quoted in Intro- 
duction, this makes the correspondence between critical 
points of the ASM and FES rather problematic. 

We consider here a nearly closed sandpile with a single 
open site io- An idea is to retrace carefully the process of 
adding particles from an initial state to the crirical state 
of the ASM. At the initial stage of the process, the prob- 
ability of dissipation via the sink at io is negligibly small 
for the large lattices and the evolution of the system al- 
most coincides with the dynamics of the FES. When the 
density approaches a critical value, the large avalanches 
appear and the probability of dissipation increases. If the 
transition to the dissipative regime is sufficiently sharp, 
we can identify the transition point with p c . Accord- 
ing to the process described in Q, the system evolves 
with further adding of particles and reaches eventually 
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the stationary density of the ASM p s . The main dif- 
ference from previous studies is in the initial conditions. 
Beside the usual initial conditions 1 < ho < 3, we con- 
sider the empty lattice ho = and show that just this 
condition leads to the transition point p c = p s . 



III. SIMULATIONS FOR ASM WITH ONE SINK 

Consider the Abelian sandpile model with one sink on a 
N x N square lattice C with periodic boundary condition 
in both directions. The initial configuration is uniform, 
hi = ho with h a = 0,1,2,3. We add particles by one 
per unit time at a random site and the resulting den- 
sity of particles is measured as a function of time. At 
the first stage of the evolution, the system loses negligi- 
bly small fraction of particles to the sink, and the den- 
sity of particles increases with time linearly. Eventually, 
the system reaches some critical density p c {ho) where it 
starts losing a macroscopic amount of particles. At this 
point, the density profile changes sharply and the den- 
sity tends monotonically to the stationary density of the 
ASM p s = 25/8. At the transition point, the density 
p(r) has a kink which allows one to determine p c {ho) 
with reasonable accuracy. 
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FIG. 1: The evolution of the dissipative ASM with one sink 
for ho — 3, number of samples - 1000, N = 500. 
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During the initial stage of evolution, the number of 
particles in the system grows linearly with time r. Rare 
dissipative events occur only if the cluster of sites with 
h = 4 is situated in a close vicinity of io and the added 
particle hits there. When the density of particles reaches 
a threshold level, a kind of percolation appears which 
can be termed a " weak percolation" . It implies that sites 
with h — 4 do not percolate yet, but additional sites 
increasing their height during an avalanche, produce to- 
gether with them a percolation set. The point of the weak 
percolation is rather low (p c 3.09 in our simulations). 
Therefore, the density of particles continues a slow in- 
crease despite the big avalanches spreading through the 
system. The weak percolation point is visible in Fig. 2 as 
a kink following the linear part of the function p(r). 



B. h = 2 

Due to lowering of the initial density, the percolation 
picture of the sites with h = 4 becomes more pronounced, 
getting close to the usual site percolation. The percola- 
tion dynamics versus the avalanche dynamics for differ- 
ent background densities in the FES has been discussed 
recently by Park [jjj • It was demonstrated that the crit- 
ical density of transition into the unstable state strongly 
depends on the background density and may differ con- 
siderably both from accepted values of p s and pth- Our 



FIG. 2: The evolution of the dissipative ASM with one sink 
for ho = 2, number of samples - 1000, N — 500. 

results for ho = 2 confirm this conclusion. The den- 
sity of particles p(r) grows strictly linearly in time up to 
the transition point p c w 3.134 which can be associated, 
due to rare dissipation events, with the critical density 
of the FES for the same initial condition. At the tran- 
sition point, the avalanche mechanism is activated and 
the dissipation reduces density towards p s . It may seem 
strange that lowering of the initial density from h = 3 to 
ho = 2 leads to the increase of the transition point from 
p c ~ 3.09 to p c 3.134. However, this is the result of the 
competition between larger number of percolation sites 
with h = 4 for ho — 2 and the reduced density of parti- 
cles outside the percolation cluster. Of course, the value 
of the transition point cannot be derived directly from 
the site-percolation critical probability P s — 0.592..., be- 
cause the distribution of heights outside the percolation 
cluster remains unknown. 



C. h = l 

This case, basically is very similar to the previous one. 
However, a tendency of increasing of the transition point 
with lowering of the initial density observed in the case 
ho = 2 is reversed. Remarkably, the value p c ~ 3.12528 
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for the relatively small lattice N — 500 is very close (af- 
ter reduction by one) to the threshold value of the FES 
p th = 2.125278 ±0.0000004 reported in for the lattice 
N = 512. This supports the idea that the nearly closed 
sandpile behaves almost identically to the FES during 
the non-dissipativc part of evolution. Indeed, the simu- 
lations in [3j were started with the initial conditions £ = 
for the sandpile model with stable occupation numbers 
C = 0, 1, 2, 3 which correspond to heights h = 1, 2,3, 4 in 
our notations. Therefore, £ = for the FES in [3j cor- 
responds to h = 1 in our simulations. Analogously, our 
ho = corresponds to £ = — 1. 



D. ho = 

The initial condition ho = is special and turns out 
to be different from the other initial conditions. The 
similarity between dynamics of the FES and the nearly 
closed ASM observed in the previous cases allows us to 
use an important statement proved for the FES: if a con- 
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FIG. 3: The evolution of the dissipative ASM with one sink 
for ho = 1, number of samples - 1000, N = 500. Two 
horizontal lines correspond to values pthiQ) = 3.125 and 
p th {l) = 3.125288.... 



3.126 



3.125 



p(T) 3.124 




3.123 



3.122 

3.122 3.13 



3.14 3.15 3.16 3.17 



3.18 



FIG. 4: The evolution of the dissipative ASM with one sink 
for ho = 0, number of samples - 1000, N = 500. Two 
horizontal lines correspond to values pth{0) = 3.125 and 
pth(l) = 3.125288.... 



figuration stabilizes after an avalanche process, there is 
some site that has not toppled [l3[ (see a similar the- 
orem in [13, HH). We may expect (and we do observe 
this in our simulations) that there is a set of sites with 
non-zero density which never topple during the evolution 
of the nearly closed ASM from the initial state up to the 
threshold value where all sites topple ultimately. Then, 
the occupation number at these sites does not affect the 
dynamics, giving at the the same time a contribution to 
the total density. If so, the threshold point p c s» 3.1252 
found for ho = 1 can be lowered for ho = to the critical 
point of the ASM. Our simulations confirm this expecta- 
tion. Adding particles to the empty lattice, we observe as 
above the linear growth of density p(r) which indicates 
almost non-dissipative character of evolution. The linear 
part changes sharply at the threshold value p c ~ 3.125 
into a function slightly fluctuating around p(r) = const. 
The coincidence of p c and p s implies that the percolation 
dynamics does not dominate anymore and the avalanche 
process approaches its standard critical behavior. 

Among the initial conditions that we have considered, 
ho = is the only one for which the threshold density 
p s appears to be close or equal to the stationary density 
25/8. We expect that this remains true for other initial 
conditions with non-positive heights. 



IV. THE THRESHOLD DENSITY OF FES 

Along with the nearly closed sandpile, we considered 
the FES for initial conditions ho = 0,1,2,3 to confirm 
the relation p c = p t h and to show the power-law conver- 
gence of the threshold density to its asymptotical value. 
To do that, we have considered N x N lattices for vari- 
ous N , and periodic boundary conditions in both direc- 
tions. On each of these lattices, we run 10 7 samples of 
fixed energy sandpile model and measure average thresh- 
old density pth(N,ho) and its standard deviation. After 
the extrapolation of pth(N, ho) we found numerically the 
asymptotic value of the threshold density and its finite- 
size corrections for large N. The accuracy 6pth(N, ho) of 
the result p t h{N, ho) is estimated via its standard devia- 
tion Dp t h(N, ho) by the following formula: 

Dpth(N) 



Sp th {N, h ) 



(13) 



V#samples 

The results of the simulations on FES with various 
initial conditions are shown in Figs. [SHH1 The sizes of 
dots in the figures are proportional to their accuracy of 
the simulations. For ho — we found 

pth(N, 0) ~ 3.1250(1) - 0.5(5)i\T* + 1.6(1).aH . (14) 
For ho = 1 we found 

p th (N, 1) ~ 3.12528(5) - 0.5(5) N~Z + l.6(l)N~i . (15) 

For ho = 2 we found 

p th (N, 2) ~ 3.13400(2) - 0.06(5)A^ + 18.5(A)N~¥ . 

(16) 
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FIG. 5: The threshold densities of the fixed energy sandpile 
for ho — 0, number of samples - 10 7 and N = 10, 11, . . . , 100. 
Two horizontal lines correspond to values pth{0) = 3.125 and 
pt h (l) =3.125288.... 




FIG. 8: The threshold densities of the fixed energy sandpile 
for hg = 3, number of samples - 10 7 and N = 10, 11, . . . , 100. 
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FIG. 6: The threshold densities of the fixed energy sandpile 
for ho — 1, number of samples - 10 7 and N = 10, 11, . . . , 100. 
Two horizontal lines correspond to values pth{0) = 3.125 and 
pth(l) =3.125288.... 
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FIG. 7: The threshold densities of the fixed energy sandpile 
for ho — 2, number of samples - 10 7 and N = 10, 11, . . . , 100. 



For ho = 3 we found 
pth(N, 3) ~ 3.0938(7) + 0.457(4)A^ 



1.010(7)iV- 



CONCLUSION 



The results of our numerical simulations support the 
following three statements: (1) the density of sand held 
by the lattice with a single dissipative site (or with dis- 
sipation rate tending to zero) , as a function of the quan- 
tity of sand deposed, shows a discontinuity at a value p c , 
which can be identified with the threshold density of the 
corresponding FES (with the same initial condition); (2) 
the critical value p c depends on the initial condition, and 
consequently, there is presumably also a range of values 
for the FES threshold density, also depending on the ini- 
tial condition and potentially on the way sand is dropped; 
(3) for the specific uniform initial configuration ho = 0, 
and probably for others as well, the critical is close or 
equal to the stationary density, p c = p s = 25/8. One 
can expect that if the system starts its evolution from 
sufficiently deep initial conditions, i.e. from initial con- 
dition with large negative heights, the stationary value 
25/8 coincides with p c and pth exactly. 



(17) 
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